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MODULATED SINE WAVES FOR 
DIFFERENTIAL ABSORPTION 
MEASUREMENTS USING A CW LASER 
SYSTEM 

5 

CROSS-REFERENCE TO RELATED PATENT 
APPLICATIONS 

This application claims the benefit of priority to and is a 
continuation-in-part of U.S. patent application Ser. No. 10 
13/167,093, filed on Jun. 23, 2011, which claims the benefit 
of priority to U.S. Provisional Patent Application Ser. No. 
61/357,687, filed on Jun. 23, 2010. This application also 
claims the benefit of priority to U.S . Provisional Patent Appli- 
cation Ser. No. 61/862,626, filed on Aug. 6, 2013. The con- 15 
tents of the foregoing applications are hereby incorporated by 
reference in their entirety. 

STATEMENT REGARDING FEDERALLY 
SPONSORED RESEARCH OR DEVELOPMENT 20 

The invention was made by employees of the United States 
Government and may be manufactured and used by or for the 
Government of the United States of America for governmen- 
tal purposes without the payment of any royalties thereon or 25 
therefore. 

FIELD OF THE INVENTION 

The present invention relates generally to range detection 30 
and differential absorption utilizing Light Detection .And 
Ranging (LiDAR), radar, sonar, or other signal. 

BACKGROUND OF THE INVENTION 

35 

The U.S. National Research Council’s report entitled Earth 
Science and Applications from Space: National Imperatives 
for the Next Decade and Beyond identified the need for a 
middle-term space mission of ASCENDS (NRC, “Earth Sci- 
ence and Applications from Space: National Imperatives for 40 
the Next Decade and Beyond,” The National Academies 
Press, Washington, D.C., 2007.) The primary objective of the 
ASCENDS mission is to make C0 2 column mixing ratio 
measurements during the day and night over all latitudes and 
all seasons and in the presence of thin or scattered clouds. 45 
These measurements would be used to significantly reduce 
the uncertainties in global estimates of CO, sources and 
sinks, to provide an increased understanding of the connec- 
tion between climate and atmospheric CO,, to improve cli- 
mate models, and to close the carbon budget for improved 50 
forecasting and policy decisions. (NRC, “Earth Science and 
Applications from Space: National Imperatives for the Next 
Decade and Beyond,” The National Academies Press, Wash- 
ington, D.C., 2007.) 

The ASCENDS mission requires an active sensor that has 55 
the needed CO, column measurement precision from space to 
determine the global variability of C0 2 in the troposphere and 
thereby determine CO, sources and sinks at the surface. The 
critical component in the ASCENDS mission is the Laser 
Absorption Spectrometer (LAS) for CO, column measure- 60 
ments with weighting towards the mid to lower troposphere 
using Integrated Path Differential Absorption (IPDA) tech- 
nique. 

IPDA lidar systems measure the difference of the total gas 
absorption along the path length between two or more laser 65 
wavelengths in the path of the laser to a target. For the earth’ s 
climate system, there are large changes in the atmospheric 


2 

variables such as pressure and temperature profiles, clouds, 
aerosols, and the profiles of water vapor and other absorbing 
gases from local spot to global scales. Also, the surface reflec- 
tance varies with many factors and variables such as surface 
type, roughness, slope, vegetation, and soil moisture. The 
variations in these atmospheric and surface factors and vari- 
ables would generate large changes in received lidar powers 
from the surface or other targets, which would introduce huge 
problems in determination of C0 2 column amounts when 
individual lidar wavelengths are investigated. The advantage 
of IPDA LAS system is that the strong effects of the variations 
in atmospheric states and surface reflection on lidar return 
powers would be eliminated from the ratio of received lidar 
powers of two close wavelengths. All the influences on the 
lidar power from these environmental variables are the same 
for the two or more close wavelengths. The ratio removes all 
environmental effects except the differential absorption opti- 
cal depth at the studied LAS wavelengths. Thus, high accu- 
racy of C0 2 column measurements (within 1 ppm) can be 
achieved from the C0 2 differential optical depth. 

Basic CO, IPDA LAS measurement approaches include 
pulsed (Grady .1. Koch, Bruce W Barnes, Mulugeta Petros, 
Jeffrey Y Beyon, Farzin Amzajerdian, Jirong Yu, Richard E 
Davis, Syed Ismail, Stephanie Vay, Michael J Kavaya, Upen- 
dra N Singh. “Coherent Differential Absorption Lidar Mea- 
surements of C02,” Applied Optics, Vol. 43 Issue 26, pp. 
5092-5099 (2004) doi: 10.1364/AO.43 .005092) and (James 
B. Abshire, Haris Riris, Graham R. Allan, Clark J. Weaver, 
Jianping Mao, Xiaoli Sun, William E. Hasselbrack, S. Ran- 
dolph Kawa “Pulsed airborne lidar measurements of atmo- 
spheric CO, column absorption,” Presented at the 8th inter- 
national carbon dioxide conference, ICDC8, in Jena 
Germany 13-19 Sep. 2009, Tellus B, 62: 770-783, doi: 
10.1 1 1 1/j. 1600 0889.2010.00502.x) andCW systems. Direct 
detection pulsed systems use a technique where measure- 
ments for each wavelength are done sequentially through 
separate columns. This works well but errors can occur if 
there are differences in reflectivity and solar background 
radiation from shot to shot and one wavelength to the next. 
Also, the atmospheric air mass measured for pulse systems is 
different from one wavelength pulse to other wavelength 
pulses. The CW approach, however, transmits lidar powers of 
different wavelengths in the same aligned laser beam, funda- 
mentally senses the same column of atmosphere and spot of 
surface, and thus, removes the potential errors for pulse sys- 
tems. Within the CW approach, IM techniques are used to 
separate the signals from different wavelengths, provide 
ranging capability, and discriminate surface returns from 
multiple other returns such as those of thin clouds, for 
instance. They are also useful for separating the laser signals 
from background solar returns, which occurs in a natural way 
through the matched filter processing techniques used. Poten- 
tial IM-CW techniques include pn code CW (Joel F. Camp- 
bell, Narasimha S. Prasad, Michael A. Flood “Pseudorandom 
noise code-based technique for thin-cloud discrimination 
with C0 2 and 0 2 absorption measurements,” Opt. Eng. 
50(12), 126002 (Nov. 18, 2011), doi:10.11 17/1.3658758)and 
(Joel F. Campbell, Michael A. Flood, Narasimha S. Prasad, 
and Wade D. Hodson “A low cost remote sensing system 
using PC and stereo equipment,” American Journal of Phys- 
ics, Vol. 79, Issue 12, pp. 1240, December 2011), linear swept 
sine wave CW (R. Agishev, B. Gross, F. Moshary, A. Giler- 
son, and S. Ahmed “Atmospheric CW-FM-LD-RR ladar for 
trace-constituent detection: a concept development,” Appl. 
Phys. B 81, 695-703(2005), doi: 10.1007/s00340-005-1919- 
x), (Oscar Batet, Federico Dios, Adolfo Comeron, and Ravil 
Agishev “Intensity-modulated linear- frequency -modulated 
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continuous-wave lidar for distributed media : 1'u ndamcn ta I s of 
technique,” Applied Optics, Vol. 49, No. 17, pp. 3369-3379, 
lOJun. 2010doi: 10.1364/A0.49.003369), (Masaharu Imaki, 
Shmnpei Kameyama, Yoshihito Hirano, Shinichi Ueno, 
Daisuke Sakaizawa, Shuji Kawakami, Masakatsu Nakajima 
“Laser absorption spectrometer using frequency chirped 
intensity modulation at 1 .57 pm wavelength for C02 mea- 
surement,” Optics Letters, Vol. 37, No. 13, pp. 2688-2690, 1 
Jul. 2012 doi: 10.1364/OL.37.002688), (Edward V. Browell, 
J. T. Dobler, S. A. Kooi, M. A. Fenn, Y. Choi, S. A. Vay, F. W. 
Harrison, B . Moore III “Airborne laser C02 column measure- 
ments: Evaluation of precision and accuracy under wide 
range of conditions,” Presented at Fall AGU Meeting, San 
Francisco, Calif., 5-9 Dec. 2011) (Edward V. Browell, J. T. 
Dobler, S. A. Kooi, M. A. Fenn, Y. Choi, S. A. Vay, F. W. 
Harrison, B. Moore III “Airborne validation of laser C02 and 
02 column measurements,” Proceedings, 1 6th Symposium on 
Meteorological Observation and Instrumentation, 92nd AMS 
Annual Meeting, New Orleans, La., 22-26 Jan. 2012. https:// 
ams.confex.com/ams/92Annual/webprogram/ 
Paperl97980.html), and (Jeremy T. Dobler, F. Wallace Har- 
rison, Edward V. Browell, Bing Lin, Doug McGregor, Susan 
Kooi, Yonghoon Choi, and Syed Ismail “Atmospheric C02 
column measurements with an airborne intensity-modulated 
continuous wave 1.57 pm fiber laser lidar,” Applied Optics, 
Vol. 52, Issue 12, pp. 2874-2892 (2013)), unswept sine wave 
CW (Songsheng Chen, Yingxin Bai, Larry B. Petway, Byron 
L. Meadows, Joel F. Campbell, Fenton W. Harrison, Edward 
V. Browell “Digital Lock-in detection for multiple-frequency 
intensity-modulated continuous wave lidar,” 26th Interna- 
tional Laser Radar Conference, SI P-38, Porto Heli, Greece, 
Porto Heli, Greece, 25-29 Jun. 2012), (Shumpei Kameyama, 
Masaharu Imaki, Yoshihito Hirano, Shinichi Ueno, Shuji 
Kawakami, Daisuke Sakaizawa, Toshiyoshi Kimura, 
Masakatsu Nakajima “Feasibility study on 1.6 mm continu- 
ous-wave modulation laser absorption spectrometer system 
for measurement of global C02 concentration from a satel- 
lite,” Applied Optics, Vol. 50, No. 14, pp. 2055-2068 (2011), 
doi: 10. 1364/AO. 50.002055) (Michael Dobbs, Jeff Pruitt, 
Nathan Blume, David Gregory, William Sharp “Matched Fil- 
ter Enhanced Fiber-based Lidar for Earth, Weather and 
Exploration,” NASA ESTO conference, June 2006 http:// 
esto.nasa.gov/ conferences/ estc2006/papers/b4p3 .pdf) . 
(Dobbs, M. E., J. Dobler, M. Braun, D. McGregor, J. Over- 
beck, B. Moore III, E. V. Browell, and T. Zaccheo “A Modu- 
lated CW Fiber Laser-Lidar Suite for the ASCENDS Mis- 
sion,” Proc. 24th International Laser Radar Conference, 
Boulder, Colo., 24-29 Jul. 2008). and (Jeremy T. Dobler, ITT 
Exelis, Fort Wayne, Ind.; and J. Nagel, V. L. Temyanko, T. S. 
Zaccheo, E. V. Browell, F. W. Harrison, and S. A. Kooi 
“Advancements in a multifunctional fiber laser lidar for mea- 
suring atmospheric C02 and 02,” Proceedings, 1 6th Sympo- 
sium on Meteorological Observation and Instrumentation, 
92nd AMS Annual Meeting New Orleans, La., 22-26 Jan. 
2012. https://ams.confex.com/ams/92Annual/webprogram/ 
Paper202790.html), non-linear swept sine wave CW (Joel F. 
Campbell “Nonlinear swept frequency technique for C02 
measurements using a CW laser system,” Applied Optics, 
Vol. 52, Issue 13, pp. 3100-3107 (2013)), and other similar 
techniques. The baseline technique used by NASA Langley 
Research Center in partnership with ITT Exelis for the 
ASCENDS mission is linear swept frequency CW (Jeremy T. 
Dobler, F. Wallace Harrison, Edward V. Browell, Bing Lin, 
Doug McGregor, Susan Kooi, Yonghoon Choi, and Syed 
Ismail “Atmospheric C02 column measurements with an 
airborne intensity-modulated continuous wave 1 .57 pm fiber 
laser lidar,” Applied Optics, Vol. 52, Issue 12, pp. 2874-2892 
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(2013)). Experiments have demonstrated the utility of this 
technique for both range and C02 column absorption mea- 
surements. An advantage of this technique is that by simulta- 
neously transmitting orthogonal online and offline signals 
5 together through a single column of atmosphere, one may 
make a simultaneous online/olfline lidar return power mea- 
surement, thereby minimizing error caused by variations in 
surface reflectivity and reflected solar background radiation. 

10 One potential issue with some modulation techniques, 
such as linear swept frequency, is the potential for side lobes, 
which could introduce small bias errors when intermediate 
scatterers are close to measurement targets. This is an advan- 
tage of digital techniques such as pseudorandom-noise (PN) 
15 codes and non-linear frequency modulation. A known tech- 
nique uses a pure maximum length (ML) sequence with time 
shifting to separate channels (Joel F. Campbell, Narasimha S. 
Prasad, Michael A. Flood “Pseudorandom noise code-based 
20 technique for thin-cloud discrimination with C02 and 02 
absorption measurements,” Opt. Eng. 50(12), 126002 (Nov. 
18, 2011), doi: 10.1 1 17/1.3658758). This technique may have 
the advantage of eliminating the side lobe issue for certain 
, s hardware configurations in LAS systems, while allowing for 
multiple channels that share the same bandwidth. Addition- 
ally, the time shifting approach makes the data processing of 
multi-channel laser returns very simple since one only needs 
to perform a single matched filter correlation for all channels 
3< 1 instead a separate matched filter correlation for each channel . 
One potential disadvantage with this technique is that not all 
hardware has a suitable low frequency response, which 
applies to the transmitter as well as the receiver. Since a ML 
35 sequence has a main power band at zero frequency, using that 
technique in such hardware will filter out the main power 
band and ruin the autocorrelation properties, which can result 
in non-orthogonality in a time shifted multi-channel configu- 
4(| ration. This problem could be solved by putting the PN code 
on a carrier, performing amplitude demodulation, and con- 
ducting matched filter correlation with the reference PN code 
to obtain the range capability. A low cost demonstration of 
this was implemented using Hilbert Transforms forms for 
45 amplitude demodulation. That system had better autocorre- 
lation properties but still exhibited minor artifacts and the 
demodulated PN code exhibited minor distortions and ring- 
ing (Joel F. Campbell, Michael A. Flood, Narasimha S. 
5Q Prasad, and Wade D. Hodson “A low cost remote sensing 
system using PC and stereo equipment,” American Journal of 
Physics, Vol. 79, Issue 12, pp. 1240, December 2011). 

BRIEF SUMMARY OF THE INVENTION 

55 

A differential absorption measurement according to the 
present invention utilizes a two-wavelength approach, which 
uses one laser with the wavelength tuned to the gas absorption 
60 line (on-line), and the other tuned away from that absorption 
line (off-line). The differential absorption optical depth for 
the CO, column may be derived from the ratio of the received 
optical power from each wavelength. 

These and other features, advantages, and objects of the 
65 present invention will be further understood and appreciated 
by those skilled in the art by reference to the following speci- 
fication, claims, and appended drawings. 



US 9,097,646 B1 


5 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. 1 is a partially fragmentary diagram of a shift register 
representation for generating m-sequence codes; 

FIG. 2 is a graph showing the result of cross correlation of 5 
code and shifted code computed with a fast F ourier transform 
for a single channel system; 

FIG. 3 is a graph of a return signal for a single channel 
system after subtracting noise; 

FIG. 4 is a graph of the return signal of FIG. 3 after corre- to 
lation; 

FIG. 5 is a graph of a return signal for a single channel 
system with lOx noise before correlation; 

FIG. 6 is a graph showing a return signal for a single 
channel system with lOx noise after correlation; 15 

FIG. 7 is a graph of a return signal utilizing two codes in a 
double ranged system with two channels before correlation; 

FIG. 8 is a graph of an example of a return signal utilizing 
two codes in a double ranged system with two channels after 
correlation; 20 

FIG. 9 is a block diagram of a CW LiDAR system with an 
AM modulated carrier; 

FIG. 10 is a block diagram of a CW LiDAR system with an 
FM modulated carrier; 

FIG. 11 is a block diagram of a CW LiDAR system that 25 
does not include a carrier; 

FIG. 12 is a baseline instrument block diagram of a LiDaR 
system according to another aspect of the present invention; 

FIG. 13 is a block diagram of a quadrature matched filter 
correlator used for demodulation of the multichannel time 30 
shifted amplitude and phase modulated sine wave modulation 
signals, wherein the quadrature matched filter correlator is 
similar to a conventional quadrature demodulator except the 
low pass filter section has been replaced by a matched filter 
correlation with the reference PN code kernel; 35 

FIG. 14A is a plot of a cosine wave modulation for an 
amplitude modulation case according to equation 15.1; 

FIG. 14B is a plot of a sine wave modulation for an ampli- 
tude modulation case according to equation 15.2; 

FIG. 15A is a plot of an autocorrelation function for a 40 
cosine wave with ML-sequence modulation with 1 cycle per 
code bit and 8 samples per code bit; 

FIG. 15B is a plot of an autocorrelation function for a sine 
wave with ML-sequence modulation with 1 cycle per code bit 
and 8 samples per code bit; 45 

FIG. 16A is a plot of a cosine wave modulation for PSK 
modulation case as in equation 1 7. 1 ; 

FIG. 16B is a plot of a sine wave modulation for PSK 
modulation case as in equation 17.2; 

FIG. 17 is a plot of the step response (signal modulation) of 50 
an Erbium-Doped Fiber Amplifier (EDFA); 

FIG. 18 is a plot of a hybrid sine wave pulse laser modu- 
lation signal, y, using the modulation wave form of equations 
17.1 and 17.2, assuming m=l; 

FIG. 19 is a plot showing an autocorrelation function for 55 
hybrid pulse modulation; 

FIG. 20 is a plot showing a comparison of frequency con- 
tent of filtered vs. unfiltered PN code; 

FIG. 21 is a plot showing a comparison between filtered PN 
code and unfiltered PN code; 60 

FIG. 22 is a plot showing an autocorrelation function for 
filtered vs. unfiltered PN code wherein the filtered case was 
correlated with an unfiltered reference; 

FIG. 23 is a plot showing filtered PN code and resulting 
amplitude modulated sine wave; 65 

FIG. 24 is a plot showing frequency spectra of filtered 
modulation vs. unfiltered modulation; 
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FIG. 25 is a plot showing a comparison of autocorrelation 
functions between a filtered modulation case and an unfil- 
tered modulation case; 

FIG. 26 is a plot showing filtered PN code and a resulting 
PSK modulated sine wave; 

FIG. 27 is a plot showing frequency spectra of filtered 
modulation vs. unfiltered modulation; 

FIG. 28 is a plot showing hybrid pulse modulation with 
filtered PN code modulation; 

FIG. 29 is a plot showing a comparison of the frequency 
spectra between filtered and unfiltered hybrid pulse case; 

FIG. 30 is a plot showing a comparison between filtered vs. 
unfiltered autocorrelation functions; 

FIG. 31 is a plot of a PN code modulated by a Jacobi theta 
function of the third kind formed by multiplying the PN code 
by the theta function. This creates a pulse for each bit of the 
PN code and is used as the transmitted laser modulation 
signal; 

FIG.32isaplot showing an example of a correlating kernel 
formed by taking twice the PN code, subtracting 1, and mul- 
tiplying the result by the Jacobi theta function of the third kind 
used in FIG. 31; and 

FIG. 33 is a plot of the resulting autocorrelation function 
for the theta function pulse-pn modulation. 

DETAILED DESCRIPTION OF PREFERRED 
EMBODIMENTS 

For purposes of description herein, terms such as “upper,” 
“lower” “right,” “left,” “rear,” “front,” “vertical,” “horizon- 
tal.” and derivatives may be utilized. However, it is to be 
understood that the invention may assume various alternative 
orientations and step sequences, except where expressly 
specified to the contrary. It is also to be understood that the 
specific devices and processes illustrated in the attached 
drawings and described in the following specification are 
simply exemplary embodiments of the inventive concepts 
defined in the appended claims. Hence, specific dimensions 
and other physical characteristics relating to the embodi- 
ments disclosed herein are not to be considered as limiting, 
unless the claims expressly state otherwise. 

In atmospheric science it is useful to use PN codes at 
multiple laser wavelengths in order to measure absorption at 
specific ranges in order to, for example, discriminate the 
return of a cloud from the ground. 

One aspect of the present invention is the use of time 
shifted PN codes for continuous wave (CW) LiDAR. In order 
to frame the requirements for use of time shifted PN codes, it 
is helpful to define the resolution and range in terms of known 
quantities. 

The resolution is given by: 


rco — — — , 

2 * bitrate 

where c is the speed of light and bitrate is the bitrate of the PN 
code. In our reference system res=3*10"8 meters/sec/ 
(2*50000/sec)=3000 meters. The range is: 

range=Af*res. (2.0) 

For instance with a PN code with a length of 256 the range 
would be 768 km, which is all the range that would be needed 
for a typical satellite-based LiDAR system. However, as dis- 
cussed in more detail below, there are advantages to making 
the code longer. 
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One known PN code is the m-sequence. The m-sequence is 
also the basis of many other PN codes. Hie m-sequences are 
advantageous because they have very good autocorrelation 
properties. M-sequences can be represented in a number of 
different ways. One known method utilizes linear feedback 5 
shift registers. An example of this implementation is shown in 
FIG. 1. In this arrangement, all additions are done using 
modulo 2 additions. 

M-sequences may also be represented using a generator 
polynomial: to 

G=gnX” + gm-l^‘+g m -2^"~ 2 + ■ ■ ■ +gP?+gl*'+go ■ (3-0) 

Here the g’s can be 1 or 0 and the sums are done using 
modulo two addition. However, g m =g 0 =l . 

These sequences are 2 m - 1 in length and only specific poly- 15 
nomials can be used for a particular length. For m=8 which 
corresponds to a code length of 2 s - 1=25 5 an allowed poly- 
nomial is: 

G=x 3 +x 6 +x fj +x 2 +l. (4.0) 20 

In order to generate the individual bits, this can be repre- 
sented by the recursion relation: 

pn„^=pn„®pn„ t 2®pn n+ 5®pn„ + 6 (5-0) 

For the initial values ofpn 0 -pn 7 any sequence of 0’ sand l’s 2:1 
can be chosen as the seed as long as they are not all 0. For a 
seedof(l, 0,1, 0,1, l,l)we generate the sequence: (1,0, 1,0, 1, 

1, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 

0 , 0 , 0 , 1 , 0 , 1 , 1 , 0 , 0 , 0 , 0 , 0 , 1 , 1 , 0 , 0 , 1 , 0 , 0 , 0 , 1 , 1 , 0 , 0 , 0 , 

0 , 1 , 1 , 0 , 1 , 1 , 1 , 1 , 1 , 1 , 0 , 1 , 1 , 1 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 1 , 30 
0 , 0 , 1 , 0 , 1 , 0 , 1 , 0 , 0 , 1 , 0 , 1 , 1 , 1 , 1 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 1 , 1 , 

0 , 0 , 1 , 1 , 0 , 0 , 0 , 1 , 1 , 0 , 1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 1 , 1 , 1 , 

0, 1, 1, 1, 1,0, 1, 1,0, 0, 1, 1, 1, 1, 1, 1, 1, 1,0, 0, 1,0, 1, 1,0. 

1,0, 1, 1,0, 1, 1, 1,0, 1,0, 1,0, 1,0, 1, 1, 1,0, 0, 1,0,0, 1, 1, 

0 , 1 , 1 , 0 , 1 , 0 , 0 , 1 , 1 , 0 , 0 , 1 , 1 , 0 , 1 , 0 , 0 , 0 , 1 , 1 , 1 , 0 , 1 , 1 , 0 , 35 
1 , 1 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 1 , 0 , 0 , 1 , 1 , 1 , 1 , 0 , 1 , 0 , 0 , 1 , 0 , 0 , 1 , 0 , 

0 , 0 , 0 , 1 , 1 , 1 , 1 , 0 , 0 , 0 , 1 , 0 , 1 , 0 , 0 , 1 , 1 , 1 , 0 , 0 , 0 , 1 , 1 , 1 , 1 ) 

The m-sequence codes have very good autocorrelation 
properties. The cross correlation between itself and a shifted 
version of itself is N when they are in sync and - 1 if they are 40 
not, where N is the length of the code. It is also understood the 
codes are AC coupled before any correlation is done by 
changing all O’ s in the code by - 1 . There are a number of ways 

of computing the correlation. One is by computing a table of: 

"45 

R{ref, data) = (6.0) 

i U 1 

— ^ ref (m) data (m + «) = DF r 1 (DFT ( ref)DFT(data )), 5Q 


where the * indicates the complex conjugate and DFT is the 
digital Fourier transform, which may be implemented as a 
fast Fourier transform (FFT). Although one may compute the 55 
sum directly, it requires something on the order of N 2 steps. 
Fast Fourier transforms are generally faster. This is because 
they only take something on the order of N log N steps. As an 
example, the case in the previous section can be taken, and the 
cross correlation between it and a code shifted 128 places to 60 
the right can be computed. FIG. 2 shows the results of cross 
correlation of code and shifted code computed with FFT. 

As a simulation, the return signal from a hard target may be 
represented as a combination of the delayed PN code in 
combination with noise. If M is the number of samples per bit 65 
of a Data Acquisition System (DAQ), and N is the length of 
the code, the strength of the signal after cross correlation will 
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be N*M when the return signal and PN code are in sync, and 
-M when they are out of sync in the absence of noise. 

As an example, 1 can be chosen as the amplitude of the 
return signal in combination with random noise 1 0 times that 
level. This example also utilizes a code that is 51 1 bits long, 
a bitrate of 50000 bits/sec, and a DAQ sample rate of 500000 
samples/sec. This gives a resolution of 3 km and a range of 
1 533 km. A code with the generating function G=x 9 +x 7 +x 5 + 
x 2 +l and a seed of (1,0, 1,0, 1,1, 1,1,1) is also chosen, as is a 
delay corresponding to 500 km. FIG. 3 shows the return 
signal for this example of a single channel system in absence 
of noise before correlation, and FIG. 4 shows the return signal 
after correlation for this example of a single channel system. 
FIG. 5 shows the same signal buried in noise that is 10 times 
the amplitude of the signal before correlation, and FIG. 6 
shows the signal of FIG. 5 after correlation. The signal can be 
readily picked out after correlation. 

The signal to noise ratio can be estimated. For purposes of 
discussion, white noise is first considered. It is possible to 
estimate the signal to noise ratio improvement before and 
after correlation. In this situation a DAQ samples the signal at: 


sample rate ( 8 . 0 ) 

M = — — 

bitrate 

samples per bit, where samplerate is the DAQ sample rate and 
bitrate is the PN code bit rate. This increases the apparent 
code length by M, which must be taken into account when 
applying the correlation calculation and the new code length 
becomes MN samples long. 

In general, white noise may be computed as a series of 
pseudorandom numbers ranging from - 1 to 1 . To generate 
pink noise, white noise can be filtered with an FFT filter using 
1/Sqrt[f] as a filter kernel, thereby giving it a power spectrum 
of 1/f. For atmospheric species detection it may be advanta- 
geous to be able to decode multiple channels. M-sequences 
(for example) may be utilized to construct a system with two 
or more channels. Since m-sequences have very good auto- 
correlation properties, shifted versions of them cross corre- 
late well except for a small CD component where they are 
uncorrelated and a large spike where they correlate. 

Also, if the LiDAR system is mounted to an artificial sat- 
ellite and the LiDAR system is directed down from space, 
there will be a final hard target (the ground) below which there 
will be no further returns possible. If it is known that there will 
be no further returns past a certain distance, the LiDAR sys- 
tem can be sized with twice the range that is otherwise 
needed, and that further range can be used to add a second 
shifted code that is shifted by half the doubled range. To get 
twice the range a PN code twice as long is used, while main- 
taining the same code bit rate. As long as there are no returns 
beyond the half range there will be no interference between 
the channels. The first half of the autocorrelation gives the 
range information for the first unshifted channel, and the 
second half of the autocorrelation gives the range information 
for the shifted second channel. The autocorrelation for both 
channels can be done in a single step utilizing, for example, 
software that has been configured to perform the necessary 
calculations. By doubling the length of the code the signal to 
noise ratio (snr) is increased. In general, the code can be as 
long as needed to add the desired number of channels. 

FIGS. 7 and 8 show an example of using two codes in a 
double ranged system with two channels before (FIG. 7) and 
after (FIG. 8) correlation. The example of FIGS. 7 and 8 
utilizes a code length of 51 1 with a sample rate of 500000 
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samples/sec and a code bit rate of 50000 bits/sec. This gives 
a total range of 51 1*3*10“5/(2*50000)=1533 km. A code 
with the generating function G=x 9 +x 7 +x 5 +x 2 +l, a seed of 
(1,0, 1,0, 1,1, 1,1,1) and a shifted version of that code shifted 
256 places to the right was used to generate the code. The 
result of this and 1 Ox noise is show in FIG. 7. This example is 
similar to the example above (FIGS. 2-6), except that two 
channels are utilized. 

In summary, two or more channels and two or more 
orthogonal PN codes are utilized, wherein each PN code is 
shifted in time. Two or more PN codes may be used on a single 
(or no) carrier. The shifting of the codes may be done such that 
there is equal spacing. The spacing between the time shifted 
sequences depends on the range. The minimum distance in 
bits must be: 


2 *br*R ( 9 . 0 ) 

n > 

c 

where n is the distance in bits, br is the bitrate, R is the 
maximum range, and c is the speed of light. Also, a code must 
be selected that is at least as long as the sum of all the 
distances. 

For instance, with two codes of length 511, the second is 
shifted to the right 256 places with respect to the first code. 
The shifting is done in such a way that as the end of the code 
is shifted past the boundary it wraps to the beginning. The 
maximum distance to the farthest target shall be 255 in that 
case. With four codes of length 511, each is shifted 1 28 places 
to the right of the preceding code. The maximum distance to 
the target in that case would be 127. In this way one may have 
multiple channels for a given code length. However, in gen- 
eral, the LiDAR/radar/sonar must be pointed at targets with a 
hard surface (such as the surface of the earth) no farther than 
the distance to the first code to avoid range-wrapping issues. 
In addition, it may be advantageous to either subtract the 
average from the code or from the data prior to processing to 
help alleviate noise/bias issues. 

With further reference to FIG. 9, a CW LiDAR system 1 
according to one aspect of the present invention utilizes an 
AM modulated carrier. The CW LiDAR system 1 includes a 
first laser 2 and a second laser 3. Although various lasers 
having different wavelengths could be utilized, in the illus- 
trated example the first and second lasers 2 and 3, respec- 
tively, may have a wavelength of about 1.57 microns. In 
general, the lasers 2 and 3 have slightly different wavelengths 
having different absorptions with respect to a particular gas or 
other substance. The first laser 2 may be selected to corre- 
spond to an absorption line of a particular chemical species. 
This depends on the chemical species to be measured. For 
example, there are numerous absorption lines at numerous 
different wavelengths. Thus, the wavelength of 1 .57 microns 
is but one example of a suitable wavelength. In general, the 
second laser 3 has a different wavelength having a relatively 
small absorption coefficient with respect to the gas to be 
measured. Furthermore, one or more additional lasers (not 
shown) having wavelengths that are not equal to the wave- 
lengths produced by the first and second lasers 2 and 3, 
respectively, may also be included in the CW LiDAR system 
1 . In general, the use of lasers having different wavelengths in 
differential absorption LiDAR is known, and the details with 
respect to selecting appropriate lasers for detecting a specific 
gas or other substance are not therefore described in detail 
herein. 
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A PN code generator 5 (see also FIG. 1) produces a first 
orthogonal PN code that is utilized to modulate the light from 
laser 2 as shown schematically at junction 6 of FIG. 9. Simi- 
larly, a PN code generator 7 modulates light from second laser 
5 3 at junction 8. The PN code generators 5 and 7 may comprise 
a device or it may comprise software. It will be understood 
that FIG. 1 may represent a device, software, or it may be used 
as a mathematical representation for generating Ml 
sequences. Light from the junctions 6 and 8 is amplified by 
to Semiconductor Optical Amplifiers (SOA’s) 10 and 11, 
respectively, to amplify the signal. The optical signal 13 from 
SOA’s 10 and 11 may be modulated to include a carrier signal 
as shown schematically in FIG. 9 by carrier 14 and junction 
12. In the example of FIG. 9, an amplitude modulated carrier 
15 is utilized. However, the inclusion of a carrier may not be 
required for all applications, and it will therefore be under- 
stood that the AM and FM modulated carriers shown in FIGS. 
9 and 10 respectively, are optional. Optical signal 15 is then 
supplied to transmit optical device 16, which produces an 
20 optical signal 18 that is directed towards a target. In FIG. 9, 
the target is shown schematically in dashed lines as a circle 
20. However, it will be understood that the target may com- 
prise clouds, gasses, the surface of a planet, or other features 
an interest. 

25 Light reflected from target 20 returns in the form of a return 

signal 22 that is gathered by a telescope 23. The light is then 
converted into an electronic fonn by a detector 24, and the 
electronic signal is then processed by a band pass filter 25 and 
AM demodulator 26. It will be understood that the telescope 
30 23. detector 24, band pass filter 25 and AM demodulator 26 
are an example of an arrangement in accordance with an 
embodiment of the present invention. However, other devices 
and configurations may also be utilized. The demodulated 
signal 28, first PN code 29, and second PN code 30 are 
35 supplied to a Digital Acquisition System (DAQ) 32, and soft- 
ware 33 provides for software cross correlation and process- 
ing according to the method/processes discussed in more 
detail above. Signal bias may be removed by software 33 as 
discussed above. 

40 With further reference to FIG. 10, an FM modulated CW 
LiDAR system 1 A according to another aspect of the present 
invention includes many components that are substantially 
the same as the AM CW LiDAR system 1 of FIG. 9. The 
components of the FM CW LiDAR system 1A of FIG. 10 that 
45 are substantially similar to corresponding components of AM 
CW LiDAR system 1 of FIG. 9 are designated with the same 
part numbers, except that the suffix “A” has been added to the 
item numbers in FIG. 10. First and second orthogonal PN 
codes are provided by devices 5A and 7 A. The orthogonal PN 
50 codes are provided to first and second Voltage Controlled 
Oscillators (VCO) 37 and 38, respectively, and light from first 
and second lasers 2A and 3 A, respectively, is modulated at the 
junctions designated 6 A and 11 A. The combined optical sig- 
nal 18A is directed towards a target 20A by transmit optics 
55 1 6A. A return signal 22 is received by a telescope 23A, and is 
then converted into an electrical signal by detector 24A. The 
signal then passes through a band pass filter 25A and an FM 
demodulator 36. The first and second orthogonal PN codes 
are provided to a DAQ 32A by devices 5A and 7A, as is the 
60 demodulated electrical return signal 28A. Software 33A then 
provides cross correlation and processing as described in 
more detail above. 

With further reference to FIG. 1 1 , a CW LiDAR system IB 
according to another aspect of the present invention is some- 
65 what similar to the AM CW LiDAR system 1 of FIG. 9 and the 
FM CW LiDAR system of FIG. 10. However, the CW LiDAR 
system IB does not utilize a carrier signal. The CW LIDAR 
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system of FIG. 11 includes first and second lasers 2B and 3B, 
and devices 5B and 7B that generate orthogonal PN codes. 
Alternatively, device 5B may generate unshifted PN code, 
and device 7B may generate PN code that is time shifted 
relative to the PN code generated by device 5B. The optical 5 
signal is amplified and transmitted towards a target 20B, and 
return signal 22B is received by a telescope 23B. A detector 
24B and high pass filter 25B provide a signal 28B to DAQ 
32B and software 33B for cross correlation and processing. 

It will be understood that FIGS. 9-11 represent examples of to 
ways to implement the present invention, but these examples 
are not intended to represent the sole ways to implement the 
invention. Various devices, software, and other such compo- 
nents may be utilized to implement the invention in LiDAR, 
radar and sonar systems. 15 

A IM-CW system 50 for C02 IPDA measurements accord- 
ing to another aspect of the present invention is illustrated in 
FIG. 12. For this system, modulated online and offline lasers 
are combined at 62, then amplified at 64 and transmitted a 
beam expander 68 as a single laser beam 72 that is reflected 20 
back at 76. A very small portion of this transmitted signal is 
split off at 66 and routed to a reference detector 82 and 
digitized simultaneously with the data scattered back from 
ground 74 and other intennediate scatterers in the atmosphere 
70 such as clouds and aerosols and received through a science 25 
detector 80 by means of a telescope 78. These reference and 
science data are further processed utilizing TIAs 84 and 86, 
respectively, optional band pass filters 88 and 90, respec- 
tively, and DAQs 92 and 94 respectively, to discriminate 
between ground and intermediate scatterers using the 30 
matched filter 96, and to obtain differential absorption power 
ratio and range estimates. The resulting data products 104 can 
be further processed, stores, and analyzed as required. 

The on-line and off-line laser modulation signals using the 
time shifted approach can be expressed as: 35 

A„y=l+mo(I), 
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due to the CO, absorption, (I is the backscatter coefficient, 6 is 
the extinction to backscatter ratio, K is a constant, where P on 
an the transmitted on-line and off-line powers, respec- 
tively. When multiple taigets including surface exist in the 
lidar path length, each taiget would generate similar received 
lidar signals like those in Equations 11.1 and 1 1 .2. The signals 
from these multiple targets are combined at the detector and 
converted to an electronic signal S(t). For an AC coupled 
receiving subsystem, we have, 

•SR) = [Ci t mcr(t -2r k !c-T) + C^mait - 2 r t /c)] ( 12 -°) 


where we sum over all k scatterers and where, 

C| S = ^P on iixp{-2e I /3(/)d/'Jexp(-2f t )exp(-2&) 


Cik = y^-PcffCX p|-2e £ /V)<f/|exp(-2&) 


where K' is a constant. These C lk and C 2k returns can be 
uniquely discriminated from other returns from different 
scatterers using the matched filter of the transmitted wave- 
form. For a ground target solve for ^ 8 ' gives: 



A oll =l+ma(t-T) (10.0) 

where m is the modulation index with a value between 40 
0 and 1 , a(t) is the repeating modulation waveform where 
-l<a(t)sl, and x is the time shifting of online modulation 
waveform related to that of offline. In general, x>2i ma Jc is 
chosen so as to avoid interference between the time shifted 
channel and intennediate scatterers, where r max is the maxi- 45 
mum distance among potential scatterers including the 
Earth’ s surface and also called unambiguous range, while c is 
the light speed. With the modulation waveforms of Equation 
10.0, the received lidar off-line and on-line optical powers 
from a single target at range r are given by 50 


Prec_off(t)= < 1U ) 

^jP^exp^-2e I" fi(/)d/'jexp(-2^)(l + mcr(r - 2r / c)) 55 


K . r , ,, (11.2) fin 

PncjnV) = — p„„ex p{-2s I fi(/)d/ ) j 

exp(— 2£)exp(— 2£)(1 + mcr(t — 2 r/c — r)) 

where £, is the total one-way column optical depth for the 65 


Where Cj and C 2g are the ground returns. Once C lg and C 2g 
are determined the column optical depth for C0 2 can be 
found. In general, this is done by cross correlating the refer- 
ence waveform with the return signal S(t), i.e., by the matched 
filter as mentioned before. This results in multiple peaks for 
these scatterers. The ground return can be identified through 
range gating (or maximum range). By choosing the reference 
modulation waveforms with perfect autocorrelation proper- 
ties (ML sequence for instance), channel separation is 
achieved though the time shifting approach since each chan- 
nel correlates at a different apparent range without interfer- 
ence between channels or scatterers. That is, each scatterer 
results in two peaks in matched filter output: one for offline, 
and the other for online with the x shifted time. The details of 
this approach are discussed elsewhere (see, e.g. Joel F. Camp- 
bell, Narasimha S. Prasad, Michael A. Flood “Pseudorandom 
noise code-based technique for thin-cloud discrimination 
with C0 2 and 0 2 absorption measurements,” Opt. Eng. 
50(12), 126002 (Nov. 18, 2011), doi: 10. 11 17/1 .3658758). 
Band Pass Filtered Modulation Waveform Generation and 
Demodulation 

The following is an improved modulation approach for 
C0 2 IPDA measurements according to one aspect of the 
present invention. This method modulates a sine wave carrier 
by an ML-sequence as in Equations 15.1 and 15.2 (below). To 
effectively represent each bit of ML-sequence in transmitted 
waveforms, and allow flexibility in the unambiguous range 
for a particular code length and sample rate, the ML-sequence 


off-line measurement that is decided by the gas absorption is oversampled by M, where M is an integer so that each code 


other than CO,, and is the one-way column optical depth bit of the ML sequence is represented by M points. Uiat bit is 



US 9,097,646 B1 


14 


13 

used to modulate a sine wave that has exactly P cycles per 
code bit. In order to satisfy the Nyquist rule we must have 
M>2P. Let pn be the original ML-sequence and PN be the 
oversampled ML-sequence such that PN(n)=pn(int((n-l)/ 
M)+l), where int is the integer part. 5 

Amplitude Modulation 

In the case of amplitude modulation, the modulation wave- 
form takes on one of two forms. 

o a (n)=PN{n)oos (2ji(k- Vz)P/M) , (15.1) 

o b {n)=PN(n)sm(2%{n-V2)P/M) (15.2) 

The 1/2 factor is introduced to make the autocorrelation func- 
tion symmetric. 

Demodulation is accomplished as shown in FIG. 13 using 15 
the reference kernel of 2PN(n)-l . This for the reference ker- 
nel is necessary to get a resulting autocorrelation function 
with off pulse values of zero. This is a quadrature demodula- 
tor (see, e.g. Joel Campbell “Synthetic quadrature phase 
detector/demodulator for Fourier transform spectrometers,” 20 
Applied Optics, Vol. 47, Issue 36, pp. 6889-6894 (2008)) with 
the low pass filter replaced by a matched filter correlation 
with the oversampled ML-sequence. The correlations may be 
efficiently computed with an FFT and take the form: 

Riref , data) = (16.0) 

j N-l 

— ^ re/ (m) data (m +n) = DFr\DFr(ref)DFT(data)), 

m= 0 

Where denotes the complex conjugate, and N is the length 

of oversampled ML-sequence. 

As an example of a modulation technique that is suitable 35 
for a band pass filtered LAS system as shown in FIG. 12, a 
special case is taken where P:=l (one cycle per code bit) and 
M=8 (8 samples per code bit). An 8th order (or 255 length) 
ML-sequence is generated using the recurrence relation pn 
(n+8)=pn(n)®pn(n+2)®pn(n+5)®pn (n+6) with the seed { 1 , 40 
0, 1, 0, 1, 1, 1, 1} for pn(l) through pn(8), where ® is the 
“exclusive or”. This is oversampled by a factor of 8. FIGS. 
14(a) and 14 (b) show a comparison between the cosine wave 
vs. sine wave modulation as in Equations 15.1 and 15.2 each 
modulation waveform produces a different autocorrelation 45 
function. Each waveform in FIGS. 14(a) and 14 (b) as well as 
the additional Figs, discussed below were plotted at 8 samples 
per code bit with linear interpolation between points. 

The autocorrelation function may be computed according 
to the block diagram of FIG. 13. The results for both the 50 
cosine and sine wave cases are shown in FIGS. 15(a) and 
15(b), respectively. The half height width is identical to the 
pure ML-sequence case. 

It can be shown that the half height width for both the 
cosine and sine wave cases is identical to the pure ML-se- 55 
quence case, though the sine wave case may be better suited 
for absorption measurements because of the more rounded 
peak. An autocorrelation function with a very sharp peak may 
have a height that is sensitive to under sampling. The sine 
wave case is less sharp than the pure ML-sequence. In both 60 
cases, and within the accuracy of numerical calculations, the 
off pulse values are approximately zero with the chosen 
parameters and normalization. The resolution of discriminat- 
ing different targets is close to the pure ML-sequence case of 
c/2Br [4,5], where Br is the code bitrate (number of code bits 65 
per second), which is related to the effective modulation 
bandwidth. 


PSK Modulation 

A type of PSK modulation can be described in the follow- 
ing way. The modulation waveform considered here takes on 
one of two forms. 

o'Jn)=(2PN(n)-l)cos(2x(n-V2)P/M), (17.1) 

& b (ny(2PN(n)-l)sm(2x(n- l A)P/M) (17.2) 

As in the previous case, the V 2 factor is introduced to make 
the autocorrelation function symmetric. This is very similar 
to the amplitude modulation case except PN(n) has been 
replaced by 2PN(n)-l . FIGS. 16(a) and 16(6) provide a com- 
parison between the cosine wave vs. sine wave modulation as 
in Equations 17.1 and 17.2, respectively. 

Demodulation is accomplished using the quadrature corr- 
elator shown in FIG. 13. This demodulation is significantly 
different than some conventional PSK demodulation tech- 
niques and is particularly useful in the time shifted channels 
approach as it allows one to separate channels without the 
need to decode the phase directly, thus allowing the use of 
multiple PSK modulated laser signals that share the same 
modulation bandwidth. According to one aspect of the 
present invention, the reference kernel PN(n) may be used. 
This reference kernel provides a resulting autocorrelation 
function with off pulse values of zero in the case of PSK. The 
autocorrelation function is identical to that shown in FIGS. 
15(a) and 15(7>) in the case where P=1 (one cycle per code bit) 
and M=8 (8 samples per code bit). In spite of identical auto- 
correlation functions, this case has a major advantage. Since 
PN(n) is used as a kernel as opposed to 2PN(n)-l , this method 
has better noise immunity. The reason is because PN(n) is 0 
for half of the length, which means we correlate over half the 
noise compared to the amplitude modulated case. As a result, 
the predicted SNR is improved by a factor of 2. 

Hybrid Pulse/CW Modulation 

.An alternative implementation utilized electronics that are 
DC coupled. This may be achieved by deleting the band pass 
filters of FIG. 12. In this case, the laser amplifier comprises an 
erbium-doped fiber amplifier (EDFA). There are at least two 
types of modulation possible with an EDFA amplifier. The 
first is pump modulation. With pump modulation there is a 
low pass filter response (see e.g. J. Freeman and J. 
Conradi “Gain modulation response of erbium-doped fiber 
amplifiers,” Photonics Technology Letters, IEEE, Vol. 5, 
Issue 2, pp. 224-226 (1993)). This would be an ideal configu- 
ration but it requires an optical amplifier for each channel. 
However, with signal modulation as has been used in the past 
(see e.g., Jeremy T. Dobler, F. Wallace Harrison, Edward V. 
Browell, Bing Lin, Doug McGregor, Susan Kooi, Yonghoon 
Choi, and Syed Ismail “Atmospheric C02 column measure- 
ments with an airborne intensity-modulated continuous wave 
1 .57 pm fiber laser lidar,” Applied Optics, Vol. 52, Issue 12, 
pp. 2874-2892 (2013)), which allows one to use the same 
amplifier for all channels, this system behaves as a lead com- 
pensator circuit ( J. Freeman and J. Conradi “Gain modulation 
response of erbium-doped fiber amplifiers,” Photonics Tech- 
nology Letters, IEEE, Vol. 5, Issue 2, pp. 224-226 (1993)), 
which is similar to a high pass filter. The step response 134 of 
such an amplifier configured for signal modulation, which is 
shown in FIG. 17, does not fully follow the input step function 
132. This is due to an excited state population depletion, 
which then settles to a steady state. 

Even if there is electronic DC coupling in the transmitter 
and receiver electronics, the nature of the optical amplifier 
used in this example to degrade the autocorrelation properties 
due to the distortion of the PN code. In order to maintain the 
integrity of the input signal the laser is continually pumped at 
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the signal level to avoid this effect. One solution is to build 
this pumping into the modulation signal. Different types of 
sine wave modulation discussed before can fulfill this pur- 
pose. An alternative is to turn the PN code on and off at a 
regular rate. This can be done with a square wave if enough 5 
bandwidth is available. A simple method is to modify Equa- 
tions 15.1 (and 15.2) to: 

cr(n) = (1 - cos(2rtn - l/2)P/ M))PN(n) - 1 = 10 

2sin 2 (rtn - 1 /2 )P/M)PN(n) - 1 


An example of this modification is shown in FIG. 18 with the 1 5 
first 100 points of the modulation wavefonn. For this 
example, P=1 and M=8 are selected. This modification has at 
least one significant advantage over the previous technique: 
for the same amplitude it uses only half the average optical 
power. The reason is because in order to transmit the previ- 
ously described modulation it must be offset. This means that 
a logical zero in that case is half the peak laser power whereas 
in this case a logical zero is zero power. 

Demodulation of this signal is relatively straight forward. 25 
Instead of demodulating the result using the block diagram in 
FIG. 13 , the return signal is correlated with 2PN-1. The 
resulting autocorrelation function is shown in FIG. 19 . 

As in the previous case, the half height width is identical to , f) 
the pure ML-sequence modulation. However, due to the more 
rounded peak of the autocorrelation function, the height mea- 
surement, which is necessary for accurate C02 absorption 
measurements, should be less sensitive to sampling than with 
the sharp peak of the pure ML-sequence case or the ML- 35 
sequence modulated sine wave case in the previous section. 

Yet the autocorrelation properties are near perfect as those in 
previous cases with the chosen parameters and normalization. 

In many LAS implementations it is desirable to limit the 
bandwidth of the PN code to reduce aliasing and thereby 411 
minimize sampling error in the measurement of the amplitude 
height of the autocorrelation function. This is provides accu- 
rate absorption measurements and makes the peak height less 
sensitive to small changes in range. However, known filters 
may not accomplish this requirement. For example, some 
known filters such as square window filter or other high order 
filter, may introduce sidelobes to the autocorrelation function. 
Gaussian filters are sometimes used in digital communica- 
tions to bandwidth limit digital data, Minimize channel bleed- 5 q 
over, and carry the maximum amount of information with the 
least bandwidth (see e.g. Conditions J. Math. Phys. 35, 1471 
(1994)). A convenient implementation is to use a periodic 
version of a Gaussian filter so that the filter can be terminated 
at the end of a cycle. Since a DFT contains the frequency 55 
content in the first half of the array, and a mirror image in the 
second half, and the PN code itself is transmitted as a repeat- 
ing waveform, a cyclic filter convolution is appropriate. 

A periodic Gaussian can be constructed in the following 
way with: 


G(x) = X exp((x - nL ) 2 / 2a 2 ) 


(19.0) 

65 


Here L is the period and a is the standard deviation. This sum 
can be done analytically as: 


^ exp((;t - nL) 2 / 2 a 2 ) = 


(20.0) 



(f 1 cos(27mx/L) 




where 0 3 is the Jacobi theta function of the third kind and the 
nome q is given by: 

q exp('-2jr J a v /t, 2 ). (21.0) 

Theta functions are widely used in the theory of elliptic func- 
tions (see e.g., John J. Benedetto, Georg Zimmermann “Sam- 
pling multipliers and the Poisson summation formula,” J. 
Fourier Anal. Appl., Vol. 3, Issue 5, pp 505-523 (1997)) and 
other applications (see e.g. Joel Campbell “The SMM model 
as a boundary value problem using the discrete diffusion 
equation,” Theor. Popul. Biol., Vol. 72, Issue 4, pp. 539-546 
(2007) and Joel Campbell “Ground state energy for the Har- 
tree-Fock equations with Dirichlet boundary”). In order to 
scale this so that is 1 at the maxima we have: 

W(x,L,q)=e 3 (3v;/L,q)IQ 1 (0,q) (22.0) 

To apply this as a filter in the frequency domain we simply 
multiply the DFT by the array W(n-1, N, q) and the inverse 
DFT is taken, where N is the length of the array and 1 <n<N. 
To show the usefulness of this filter we first show how to limit 
the bandwidth of an MLsequence and show how the autocor- 
relation properties and resolution are preserved. That band- 
width limited sequence is then used in the modulation 
examples discussed in the previous sections to show both the 
effect of bandwidth limiting and preservation of the near 
perfect autocorrelation properties for each case with only a 
minor degradation in resolution. 

Theta Function Filtered ML-Sequence 

First, take the simple case of a pure ML sequence with 
length 255 and oversampled by a factor of 8 as in the previous 
section. This is filtered with the theta function filter with 
q=0.9 and N=255x8=2040. A plot of the unfiltered vs. filtered 
PN code frequency content is shown in FIG. 20 . Filtered pn 
code bandwidth is well under the Nyquist sampling rate dem- 
onstrating aliasing is no longer an issue. 

A plot of unfiltered vs. filter PN code is shown in FIG. 21. 
This is the result of a circular convolution between the cyclic 
theta function filter and the cyclic PN code. 

FIG. 22 is the result of the correlation between 2PN- 1 and 
both the filtered and unfiltered PN code. In spite of the strong 
filtering, the pulse width is very similar, and the peak is more 
rounded, which is a desirable trait. The off pulse values are 
zero to within numerical precision. 

Amplitude Modulations Case Using Theta Function Filtered 
PN Code 

This case uses the filter PN code described earlier to modu- 
late a sine wave carrier as in a b of Equations 15.1 and 15.2 
with P=2 (2 cycles per code bit) and M=8 (8 samples per code 
bit). A plot of the resulting modulation is shown in FIG. 23. 

FIG. 24 is a lot of the frequency spectra of the modulated 
sine wave. This shows the resulting frequency content of the 
filtered modulation is bounded within the Nyquist range with 
no aliasing. 

FIG. 25 is a plot of the autocorrelation function of the 
filtered vs. unfiltered pn modulation. This shows the modu- 
lation produces more robust matched filter output for peak 
determination by producing a more rounded peak, which is 
desirable for amplitude measurements in differential absorp- 
tion estimation because it produces less sampling error. The 
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resolution degradation is very minor, and off pulse values are 
zero to within numerical precision. 

PSK Modulation Case Using Theta Function Filtered PN 
Code 

5 

This case uses the filtered PN code described earlier to 
modulate a sine wave carrier as in a' b of Equation 18.0 with 
P=2 (2 cycles per code bit) and M=8 (8 samples per code bit). 

A plot of the resulting modulation is shown in FIG. 26. 

The plot of the autocorrelation function of the filtered vs 10 
unfiltered pn PSK modulation is identical to FIG. 25 of the 
previous section. 

Hybrid Pulse Modulation with Theta Fimction Filtered ML 
Sequence 15 

We take the case using Equations 17.1 and 17.2 with P=1 
(one cycle per code bit) and M=8 (8 samples per code bit), and 
use the same filtered PN code as in the previous two examples 
with q=9.0 FIG. 28 shows the resulting modulation. 

20 

Hie Resulting frequency spectra are shown in FIG. 29. The 
original modulation with the chosen parameters already does 
a good job at limiting the high frequencies as one might 
deduce from FIG. 17. This is because there are no sharp 
transitions as there are with the pure ML-sequence. The fil- 25 
tered version limits the spectra even more. Nyquest rate, but 
the unfiltered modulation has a limited bandwidth. 

FIG. 30 is a plot of the autocorrelation function of the 
filtered vs unfiltered pn code modulation. Both cases show 
very good autocorrelation properties, and off pulse values are 11 1 
zero to within numerical precision. 

The three different PN modulation techniques described 
above can be used with different types of hardware. These 
include the baseline system (FIG. 12) witha signal modulated 35 
EDFA, the baseline system (FIG. 12) with the bandpass filters 
deleted, or the baseline system (FIG. 12) with the bandpass 
filters deleted and using separate EDFAs with pump modu- 
lation. In the case of DC coupled electronics combined with a 
laser amplifier with pump modulation, a pure ML sequence is 40 
ideal. In the case of the baseline, although the amplihide 
modulated sine wave version works well, the PSK modula- 
tion case discussed above may be superior in that it has 
superior signal to noise. If the electronics is DC coupled but 
the laser has a transient response such as the case with a signal 4> 
modulated EDFA, the hybrid pulse version appears to be best. 

The new cyclic filter based on the Jacobi theta function 
discussed above is useful for bandwidth limiting digital pn 
codes without seriously degrading the resolution and result- 50 
ing in autocorrelation properties completely without side- 
lobes or other artifacts. A primary benefit is the reduction of 
aliasing, which reduces sampling error when measuring the 
amplitude of the autocorrelation fimction. 

Appendix A: Periodic Gaussian as a Jacobi Theta Function 55 

The sum represented by Equation 19.0 may be performed 
analytically using the Poisson summation formula “(see e.g. 
John J. Benedetto, Georg Zimmermann Sampling multipliers 
and the Poisson su mm ation formula J. Fourier Anal. Appl., 6Q 
Vol. 3, Issue 5, pp 505-523 (1997))”, which states: 


V m = V m 


where 


»-£ 


/ («)exp( - Inikn) dn 


(A2) 


is the Fourier Transform of f. In the case of Equation 1 0.0 this 
becomes: 


m=£ 


exp((;t - nL) 2 / 2cr 2 )exp(—2mkn) dn — 


(A3) 


V In | — |exp(— 2x 2 cr 2 k 2 / L 2 )QXp(-2nkx / L) 


so that 


F{x) = V 2tt j — | ^ exp(-27r 2 cr 2 k 2 / L 2 )exp(-2nkx / L) = 

k=- 00 

y/2jr | — | 1 + 2 U cos(2^/.'.y/ L) = \!2^\jjd 3 !nxlL. q) 


(A4) 


Another aspect of the present invention is a Gaussian pulse 
PN LiDAR for differential absorption measurements. The 
Gaussian pulse PN LiDAR may utilize a LiDAR system 50 as 
shown in the blocked diagram of FIG. 12 except that the 
Gaussian pulse PN LiDAR does not require band pass filters 
88 and 90. The Gaussian pulse PN LiDAR may otherwise be 
substantially identical to the system 50 of FIG. 12. 

The Gaussian pulse LiDAR system transmits a periodic 
Gaussian pulse that is turned on and off in a very specific way 
based on a digital sequence with good auto correlation prop- 
erties such as an ML-sequence, Barker code, etc. that have an 
amplitude from 0 to 1 . The code is oversampled so that one bit 
extends over M samples. Each code bit is replaced by a 
Gaussian pulse, and this is what is transmitted through the 
LiDAR. A correlating kernel for that code is generated by 
taking that code, multiplying it by 2 and subtracting 1 , then 
oversampling by a factor of q as before. Each bit of the 
correlating kernel is also replaced by Gaussian pulse that may 
or may not have the same width as the transmitted Gaussian 
pulses. An approximation of this system canbe constructed as 
follows. Let PN be the digital sequence (some type of PN 
code) and PN be the oversampled PN code. The oversampled 
code is then: 

PN(n)=pn(int((n-l/M)+l ) (23.0) 

where int is the integer part. We define a fimction W such 
that 

W(x,q)=Q 3 (nx,q)/Q 3 (0,q) (24.0) 

where 0 3 is the Jacobi theta function of the third kind and q is 
the nome. By definition this function is derived as a sum of 
Gaussian pulses and q is related to the standard deviation of 
these pulse by, 

q exp(— 2jr 2 <F) (25.0) 

where a is the standard deviation of the Gaussian pulse. 

To generate the outgoing pulses we write 

Sou,(n)=PN(n ) W((n-'k)!M-V 2 ,q) (26.0) 

The correlating kernel is generated from 

K(ri)-(2PN{n )- 1 ) W( (n-Vi)IM-V 2 ,q ' ) (27.0) 

Where q' may or may not be the same as q. Choosing a 
different (usually smaller) value for q prime is useful in 
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adjusting the shape of the autocorrelation function to make it 
more suitable for absorption measurements. 

An example plot of the outgoing signal is shown in FIG. 31. 
An example of the correlating kernel is shown in FIG. 32. 
The resulting cross correlation between the outgoing signal 
and correlating kernel is shown in FIG. 33. 

It is to be understood that variations and modifications can 
be made on the aforementioned structure without departing 
from the concepts of the present invention, and further it is to 
be understood that such concepts are intended to be covered 
by the following claims unless these claims by their language 
expressly state otherwise. 

What is claimed is: 

1. A method of utilizing an intensity modulated laser beam 
to provide for measurement of differential absorption in a 
continuous wave lidar system, the method comprising: 

utilizing a laser to provide a first light beam having a first 
wavelength; 

utilizing a laser to provide a second light beam having a 
second wavelength that is not equal to the first wave- 
length; 

utilizing a PN code to provide a modulation signal having 
an intensity modulated (IM) repeating wave fonn corre- 
sponding to the PN code; 

providing intensity modulation of the first and second light 
beams whereby the first and second light beams com- 
prise first and second transmitted signals corresponding 
to the modulation signal, wherein the first transmitted 
signal is time shifted relative to the second transmitted 
signal; 

directing the first and second light beams at a target such 
that at least a portion of the first light beam is reflected, 
and at least a portion of the second light beam is 
reflected; 

receiving at least a portion of the reflected first and second 
light beams to form a received signal; 

demodulating the received signal utilizing the modulation 
signal to provide first and second distinct amplitude 
peaks; 

determining a magnitude of the first amplitude peak; 

determining a magnitude of the second amplitude peak; 
and 

utilizing the magnitudes of the first and second peaks to 
determine a relative absorption of the first and second 
light beams. 

2. The method of claim 1, wherein: 

the IM repeating wave fonn comprises an amplitude modu- 
lated sine wave. 

3. The method of claim 2, wherein: 

the received signal is demodulated utilizing a quadrature 
matched filter correlator. 

4. The method of claim 3, wherein: 

the quadrature matched filter correlator provides a matched 
filter correlation utilizing a PN code kernel. 

5. The method of claim 1, wherein: 

the IM repeating waveform comprises a sine wave ampli- 
tude modulated by a PN code. 

6. The method of claim 5, wherein: 

the received signal is demodulated utilizing a quadrature 
matched filter correlator. 

7. The method of claim 1, wherein: 

the IM repeating wave fonn comprises a phase modulated 
sine wave. 

8. The method of claim 7, wherein: 

the received signal is demodulated utilizing a quadrature 
matched filter correlator. 
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9. The method of claim 1, wherein: 

the IM repeating wavefonn comprises a PN code phase 
modulated sine wave. 

10. The method of claim 9, wherein: 

5 the received signal is demodulated utilizing a quadrature 
matched filter correlator. 

11. The method of claim 1, wherein: 

the IM repeating waveform comprises a frequency modu- 
lated sine wave. 

111 12. The method of claim 11, wherein: 

the received signal is demodulated utilizing a quadrature 
matched filter correlator. 

13. The method of claim 1, wherein: 

15 the IM repeating wavefonn comprises a PN code fre- 
quency modulated sine wave. 

14. The method of claim 13, wherein: 

the received signal is demodulated utilizing a quadrature 
matched filter correlator. 

20 15. The method of claim 1, wherein: 

the IM repeating waveform comprises a sine squared 
modulated PN code fonning a single pulse for each bit of 
the PN code. 

16. The method of claim 15, wherein: 

25 the IM repeating wavefonn is demodulated using a corre- 
lating kernel fonned from twice the PN code minus one. 

17. The method of claim 15, wherein: 

the IM repeating waveform is demodulated using a corre- 
lating kernel fonned from twice the PN code minus one 

30 multiplied by a sine squared function forming a single 
pulse for each bit of the PN code. 

18. The method of claim 1, wherein: 

the IM repeating waveform comprises a repeating Gauss- 
ian pulse modulated PN code fonning a single pulse for 

35 each bit of the PN code. 

19. The method of claim 18, wherein: 

the IM repeating waveform is demodulated using a corre- 
lating kernel fonned from twice the PN code minus one. 

20. The method of claim 18, wherein: 

40 the IM repeating waveform is demodulated using a corre- 
lating kernel formed from twice the PN code minus one 
multiplied by a repeating Gaussian pulse forming a 
single pulse for each bit of the PN code. 

21. The method of claim 1, wherein: 

45 the IM repeating wavefonn comprises a Jacobi theta func- 
tion of the third kind modulated PN code forming a 
single pulse for each bit of the PN code. 

22. The method of claim 21, wherein: 

the IM repeating wavefonn is demodulated using a corre- 

50 lating kernel fonned from twice the PN code minus one. 

23. The method of claim 21, wherein: 

the IM repeating wavefonn is demodulated using a corre- 
lating kernel fonned from twice the PN code minus one 
multiplied by a Jacobi function of the third forming a 

55 single pulse for each bit of the PN code. 

24. The method of claim 21, wherein: 

the IM repeating wavefonn is demodulated using a corre- 
lating kernel fonned from twice the PN code minus one 
multiplied by a repeating analog pulse fonning a single 

60 pulse for each bit of the PN code. 

25. The method of claim 1, wherein: 

the IM repeating waveform comprises an analog pulse 
modulated PN code fonning a single pulse for each bit of 
the PN code. 

65 26. The method of claim 25, wherein: 

the IM repeating wavefonn is demodulated using a corre- 
lating kernel fonned from twice the PN code minus one. 
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27. The method of claim 1, wherein: 

the PN code is filtered utilizing a cyclic frequency attenu- 
ation function based on a Jacobi theta function of the 
third kind. 

28. The method of claim 1, wherein: 5 

the PN code comprises a filtered PN code. 



